######################
#  Replication code for 'Mediating the Electoral Connection', forthcoming in the JOP
#  John Henderson and John Brooks
#  12/7/2015    
######################    

# figureVI.R
#  :: runs simulation analysis looking at attenuation given the imputation of 484 losing incumbents using 
#    their previous ideal points

rm(list=ls())
setwd('~/Dropbox/rainReplication')

fes.type=3 
non.missings = 4

source('prelimRobust.R')       

# FIGURES; 

vec=c(0.005+0.025*0:16)

modr_mat_1=matrix(NA,17,3)
modr_mat_2=matrix(NA,17,3)  
extr_mat_1=matrix(NA,17,3)
extr_mat_2=matrix(NA,17,3)
   
main_fe_dose=main_fe_vote=array(NA,nrow(rain_data))  
vts=lm(vote~as.factor(fe_id_num),data=covs,subset=full)$res
dse=lm(dose~as.factor(fe_id_num),data=covs,subset=full)$res

main_fe_dose[as.numeric(names(dse))]=dse
main_fe_vote[as.numeric(names(vts))]=vts          

covs$main_fe_dose=main_fe_dose
covs$main_fe_vote=main_fe_vote

year=rain_data$year         
lower=rain_data$lower

for(i in 1:length(vec)){
	modr=(year>1954 & lower==1 & (dose_prv>= -vec[i] & dose_prv<= vec[i])) #**
	extr=(year>1954 & lower==1 & (dose_prv<= -vec[i] | dose_prv>= vec[i])) #**	  
	
	modr=which(modr)
	extr=which(extr)
	
	mmain_iv1_fe=ivreg(main_fe_vote~
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	main_fe_dose + 
	dose_prv + vote_prv,
	~ dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	dose_prv + vote_prv + rain_day+rain_day_prev,
	subset=modr,data=covs) 
  
	mmain_iv2_fe=ivreg(main_fe_vote~
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	main_fe_dose + 
	dose_prv + vote_prv,
	~ dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	dose_prv + vote_prv + rain_weekend+rain_weekend_prev,
	subset=modr,data=covs)

	mmain_iv1_fe_sumcl=coeftest.cluster(covs[modr,],mmain_iv1_fe,cluster1='as.factor(fe_id_num)')
	mmain_iv2_fe_sumcl=coeftest.cluster(covs[modr,],mmain_iv2_fe,cluster1='as.factor(fe_id_num)')
	mmain_iv1_fe_sum=summary(mmain_iv1_fe)
	mmain_iv2_fe_sum=summary(mmain_iv2_fe)
	
	emain_iv1_fe=ivreg(main_fe_vote~
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	main_fe_dose +
	dose_prv + vote_prv,
	~ dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq +
	dose_prv + vote_prv + rain_day+rain_day_prev,
	subset=extr,data=covs) 

	emain_iv2_fe=ivreg(main_fe_vote~
	dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	main_fe_dose +
	dose_prv + vote_prv,
	~ dist_prev + midterm + pres_party + 
	black + construction + educ + 
	minc + farmer + forborn + gvtwkr + manuf + pop + unempld + 
	urban + retail + sos + gov + comp_cq + 
	dose_prv + vote_prv + rain_weekend+rain_weekend_prev,
	subset=extr,data=covs)

	emain_iv1_fe_sumcl=coeftest.cluster(covs[extr,],emain_iv1_fe,cluster1='as.factor(fe_id_num)')
	emain_iv2_fe_sumcl=coeftest.cluster(covs[extr,],emain_iv2_fe,cluster1='as.factor(fe_id_num)')
	emain_iv1_fe_sum=summary(emain_iv1_fe)
	emain_iv2_fe_sum=summary(emain_iv2_fe)

	modr_mat_1[i,1]=mmain_iv1_fe_sumcl[which(names(mmain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),1]#/
	modr_mat_1[i,2]=mmain_iv1_fe_sumcl[which(names(mmain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),2]	
	modr_mat_1[i,3]=mmain_iv1_fe_sumcl[which(names(mmain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),4]
	
	modr_mat_2[i,1]=mmain_iv2_fe_sumcl[which(names(mmain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),1]#/
	modr_mat_2[i,2]=mmain_iv2_fe_sumcl[which(names(mmain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),2]	
	modr_mat_2[i,3]=mmain_iv2_fe_sumcl[which(names(mmain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),4]	
  	 
	extr_mat_1[i,1]=emain_iv1_fe_sumcl[which(names(emain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),1]#/
	extr_mat_1[i,2]=emain_iv1_fe_sumcl[which(names(emain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),2]	
	extr_mat_1[i,3]=emain_iv1_fe_sumcl[which(names(emain_iv1_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),4]
	
	extr_mat_2[i,1]=emain_iv2_fe_sumcl[which(names(emain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),1]#/
	extr_mat_2[i,2]=emain_iv2_fe_sumcl[which(names(emain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),2]	
	extr_mat_2[i,3]=emain_iv2_fe_sumcl[which(names(emain_iv2_fe_sumcl[,1])=="dose" | names(mmain_iv1_fe_sumcl[,1])=="main_fe_dose"),4]
	
}             
 
extr_f_1=extr_f_2=modr_f_1=modr_f_2=matrix(NA,17,1000)    
v=m1=m2=e1=e2=c()   
set.seed(1005)
for(i in 1:length(vec)){            
	modr_f_1[i,]=rnorm(mean=modr_mat_1[i,1],sd=modr_mat_1[i,2],n=1000)
	modr_f_2[i,]=rnorm(mean=modr_mat_2[i,1],sd=modr_mat_2[i,2],n=1000)
	extr_f_1[i,]=rnorm(mean=extr_mat_1[i,1],sd=extr_mat_1[i,2],n=1000)
	extr_f_2[i,]=rnorm(mean=extr_mat_2[i,1],sd=extr_mat_2[i,2],n=1000)	
	v=c(v,rep(vec[i],1000))
	m1=c(m1,modr_f_1[i,])
	m2=c(m2,modr_f_2[i,])	
	e1=c(e1,extr_f_1[i,])
	e2=c(e2,extr_f_2[i,])	
}


pdf(file='simulations/simruns/comp_slide_extr.pdf')
oldpar <- par(mar = c(5.1, 4.1, 4.1, 4.1))  
hist(axes=F,breaks=200,-abs(covs$dose[which(!is.na(covs$dose)&rain_data$d_inc==1&rain_data$d_win==0 & year>1954&lower==1 | rain_data$r_inc==1&rain_data$r_win==0 & year>1954&lower==1 & !is.na(covs$dose))]),
	ylim=c(-20,20),xlim=c(-.4,0),xlab='Previous Vote Margin (Absolute)',ylab='Instrumental Estimate',main='')
j=3;k=2 
    axis(4,)
    axis(1,at=c(-.4,-.3,-.2,-.1,0),label=as.character(c(0.4,0.3,0.2,0.1,0.0))) 
    axis(2,at=c(-20,-10,0,10,20),label=as.character(c(-20,-10,0,10,20)/5))   
	mtext("Frequency of Imputed by Defeat", side = 4, line = 3)  
x1=(lowess(f=1/k,tapply(v,v,mean),j*tapply(e1,v,quantile,prob=.05)))$x
x2=(lowess(f=1/k,tapply(v,v,mean),j*tapply(e1,v,quantile,prob=.95)))$x
y1=(lowess(f=1/k,tapply(v,v,mean),j*tapply(e1,v,quantile,prob=.05)))$y
y2=(lowess(f=1/k,tapply(v,v,mean),j*tapply(e1,v,quantile,prob=.95)))$y
polygon(x=-c(rev(x1),(x1)),y=c(rev(y1),(y2)),density=35,border=F,col='slategrey')       
dev.off()       
     
pdf(file='simulations/simruns/comp_slide_modr.pdf')
oldpar <- par(mar = c(5.1, 4.1, 4.1, 4.1))  
hist(axes=F,breaks=200,-abs(covs$dose[which(!is.na(covs$dose)&rain_data$d_inc==1&rain_data$d_win==0 & year>1954&lower==1 | rain_data$r_inc==1&rain_data$r_win==0 & year>1954&lower==1 & !is.na(covs$dose))]),
	ylim=c(-20,20),xlim=c(-.4,0),xlab='Previous Vote Margin (Absolute)',ylab='Instrumental Estimate',main='')
j=3;k=2 
    axis(4,)
    axis(1,at=c(-.4,-.3,-.2,-.1,0),label=as.character(c(0.4,0.3,0.2,0.1,0.0))) 
    axis(2,at=c(-20,-10,0,10,20),label=as.character(c(-20,-10,0,10,20)/5))   
	mtext("Frequency of Imputed by Defeat", side = 4, line = 3)
x1=(lowess(f=1/k,tapply(v,v,mean),j*tapply(m1,v,quantile,prob=.05)))$x
x2=(lowess(f=1/k,tapply(v,v,mean),j*tapply(m1,v,quantile,prob=.95)))$x
y1=(lowess(f=1/k,tapply(v,v,mean),j*tapply(m1,v,quantile,prob=.05)))$y
y2=(lowess(f=1/k,tapply(v,v,mean),j*tapply(m1,v,quantile,prob=.95)))$y                                                                        
polygon(x=-c(rev(x1),(x1)),y=c(rev(y1),(y2)),density=35,border=F,col='slategrey')          
dev.off()

# end     